suppressPackageStartupMessages({
library("EnhancedVolcano")
library("limma")
library("biomaRt")
library("annotables")
library("tximport")
library("apeglm")
library("eulerr")
library("DESeq2")
library("HGNChelper")
library("tictoc")
library("DESeq2")
library("kableExtra")
library("beeswarm")
library("missMethyl")
library("gridExtra")
library("png")
library("metafor")
library("ggplot2")
library("purrr")
library("metafor")
  library("AnnotationDbi")

library("readxl")
library("ggplot2")
library("tidyverse")
library("magrittr")
library("readr")
library("eulerr")
library("RColorBrewer")
library("e1071")
library("plotly") 
library("pheatmap")
library("msigdbr")
library("mitch")
  library("rtracklayer")
  library("dplyr")
  library("org.Mm.eg.db")
  library("fgsea")
  
})

CORES=16

Read the sample manifest and alignment summary

I ran two types of aligment with kallisto (with and without ncRNA), I have comapred the result of the alingement and choosed the cDNA+ncRNA as the aligment since the aligment percentage is higher in this combo

base<- "/mnt/vol1/Mouse_model_RNA_Seq/kallisto_results_plus_ncRNA/"
idx<- "/mnt/vol1/Mouse_model_RNA_Seq/index"

samples<- read_tsv(file.path(base,"samples.txt"),col_names = "sample", show_col_types = FALSE)
samples$path <- file.path(base, samples$sample, "abundance.tsv")

align<- read_delim("/mnt/vol1/Mouse_model_RNA_Seq/kallisto_results_plus_ncRNA/alignment_summary.tsv")

old_align<-read_csv("/mnt/vol1/Mouse_model_RNA_Seq/kallisto_results/aligment_no_ncrna.csv", show_col_types = FALSE)

new_df<- align%>% mutate(index = "cDNA + ncRNA", rate = 100 * pseudoaligned/processed) %>%
 select(sample, index, rate)
old_df <- old_align %>% mutate(index = "cDNA only",                               rate = 100 * pseudoaligned/processed) %>%
  select(sample, index, rate)

df<- bind_rows(new_df, old_df)


delta<- df %>% pivot_wider(names_from = index, values_from = rate) %>%
  mutate(delta = `cDNA + ncRNA` - `cDNA only`) %>%
  arrange(delta)
df$sample <- factor(df$sample, levels = delta$sample)

ggplot(df, aes(sample, rate, fill = index)) +
  geom_col(position = position_dodge(width = 0.75), width = 0.7) +
  labs(title = "Pseudoalignment rate by sample",
       x = "Sample", y = "Rate (%)", fill = "Run") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))

Import transcripts to gene level via trimport package

tx2gene<- read.csv(file.path("/mnt/vol1/Mouse_model_RNA_Seq/tx2gene_r115.csv"), stringsAsFactors = FALSE)
stopifnot(all(c("TXNAME","GENEID") %in% names(tx2gene)))

files<- file.path(base, samples$sample, "abundance.tsv")
names(files) <- samples$sample
stopifnot(all(file.exists(files)))

txi<- tximport(files,
                type = "kallisto",
                tx2gene = tx2gene[, c("TXNAME","GENEID")],
                countsFromAbundance = "lengthScaledTPM",
                ignoreTxVersion = TRUE)

txi$counts[1:2, 1:6]

##                       BB329    BB330    BB331    BB332    BB334    BB335
## ENSMUSG00000000001 3423.529 4452.743 3423.237 3360.252 4088.225 2816.038
## ENSMUSG00000000003    0.000    0.000    0.000    0.000    0.000    0.000

Reading meta data that comes with the samples

pheno<- read_excel("/mnt/vol1/Mouse_model_RNA_Seq/minipump.xlsx",skip =1 )
colnames(pheno)

## [1] "Sample"                     "Sex"                       
## [3] "Intervention"               "RNA concentration\r\nng/ul"
## [5] "RIN"                        "Well Number"

coldata<- pheno %>%
  filter(Sample %in% samples$sample) %>%
  arrange(match(Sample, samples$sample)) %>%
  mutate(
    Sex = factor(Sex, levels = c("F","M")), # F as baseline
    Treatment = factor(Intervention, levels = c("Saline","Ang-II"))) # saline as baseline

rownames(coldata)<- coldata$Sample

Making DESeqDataSet from tximport

dds<- DESeqDataSetFromTximport(txi, colData = coldata, design = ~ Sex + Treatment + Sex:Treatment)

Filtering for low count and normalisation of the counts

keep<- rowSums(counts(dds)) >= 10
dds<- dds[keep, ]

dds<- estimateSizeFactors(dds)
norm_counts<- counts(dds, normalized = TRUE)
glimpse(norm_counts)

##  num [1:41594, 1:17] 3210.9 157.5 1444.1 196.9 10.3 ...
##  - attr(*, "dimnames")=List of 2
##   ..$ : chr [1:41594] "ENSMUSG00000000001" "ENSMUSG00000000028" "ENSMUSG00000000031" "ENSMUSG00000000037" ...
##   ..$ : chr [1:17] "BB329" "BB330" "BB331" "BB332" ...

Unsupervised clustering analyses

vsd<- vst(dds, blind = TRUE)

ann<- as.data.frame(colData(dds))[, c("Sex","Treatment"), drop = FALSE]

#sample–sample correlation heatmap
pheatmap(cor(assay(vsd)),
         annotation_col = ann,
         annotation_row = ann)

# PCA (colour = Treatment, shape = sex)
plotPCA(vsd, intgroup = c("Treatment","Sex"))

## Fit model (set baselines)
dds$Sex<- relevel(factor(dds$Sex), "F")
dds$Treatment<- relevel(factor(dds$Treatment), "Saline")
design(dds)<- ~ Sex + Treatment + Sex:Treatment
dds<- DESeq(dds)

resultsNames(dds)

## [1] "Intercept"                  "Sex_M_vs_F"                
## [3] "Treatment_Ang.II_vs_Saline" "SexM.TreatmentAng.II"

# baseline (saline): M vs F
coef_sex<- "Sex_M_vs_F"   

# Ang-II vs Saline in females (F is baseline sex)
coef_trt<- "Treatment_Ang.II_vs_Saline"  

# (Ang-II effect in M) & (Ang-II effect in F)
coef_int<- "SexM.TreatmentAng.II"

RAW results vs Shrinking the log2 fold changes

With LFC filtering

### RAW results ###

# 1. Within each Sex
# 1.1 Female: Ang-II vs Saline  
res_F_raw<- results(dds, name = coef_trt)
# 1.2 Male: Ang-II vs Saline  
res_M_raw<- results(dds, contrast = list(c(coef_trt, coef_int)))

# 2. Effect of Ang-II treatment between sex
res_int_raw<- results(dds, name = coef_int)

# 3. Baseline (saline): M vs F
res_bsl_raw<- results(dds, name = coef_sex)

# Shrinkage results
res_F_shr<- lfcShrink(dds, coef = coef_trt, type = "ashr")
res_M_shr<- lfcShrink(dds, contrast= list(c(coef_trt, coef_int)), type = "ashr")
res_int_shr<- lfcShrink(dds, coef = coef_int, type = "ashr")
res_bsl_shr<- lfcShrink(dds, coef = coef_sex, type = "ashr")

make_df<- function(raw, shr, label){
  data.frame(
    gene = rownames(raw),
    lfc_raw= raw$log2FoldChange,
    lfc_shr = shr$log2FoldChange,
    FDR_raw = raw$padj,           
    contrast = label,
    stringsAsFactors = FALSE)
}

res_all<- bind_rows(
  make_df(res_F_raw, res_F_shr, "Female: Ang-II vs Saline"),
  make_df(res_M_raw, res_M_shr, "Male: Ang-II vs Saline"),
  make_df(res_int_raw, res_int_shr,"Interaction: (M & F) treatment effect"),
  make_df(res_bsl_raw, res_bsl_shr, "Baseline (Saline): M vs F")
) 

# Raw vs Shrunken
res_long<- res_all %>%
  mutate(
    hit_raw= !is.na(FDR_raw) & FDR_raw<= 0.05 & abs(lfc_raw)>= 0.5,
    hit_shr= !is.na(FDR_raw) & FDR_raw<= 0.05 & abs(lfc_shr)>= 0.5
  ) %>%
  dplyr::select(gene, contrast, FDR_raw, lfc_raw, lfc_shr, hit_raw, hit_shr) %>%
  tidyr::pivot_longer(
    cols = c(lfc_raw, lfc_shr, hit_raw, hit_shr),
    names_to= c(".value", "lfc_type"),
    names_sep= "_"              
  ) %>%
  mutate(
    lfc_type = dplyr::recode(lfc_type,
      raw = "Raw LFC",
      shr = "Shrunken LFC"),
    lfc_type = factor(lfc_type, levels = c("Raw LFC", "Shrunken LFC")),
    contrast = factor(
      contrast,
      levels = c(
        "Female: Ang-II vs Saline",
        "Male: Ang-II vs Saline",
        "Interaction: (M & F) treatment effect",
        "Baseline (Saline): M vs F"
      )
    )
  )


x_max<- quantile(abs(res_long$lfc), 0.995, na.rm = TRUE)
x_lim<- c(-max(1.5, x_max), max(1.5, x_max))
y_max<- quantile(-log10(res_long$FDR_raw), 0.995, na.rm = TRUE)
y_lim<- c(0, max(5, y_max))

ggplot(res_long, aes(lfc, -log10(FDR_raw), colour = hit)) +
  geom_point(alpha = 0.6, size = 1.2) +
  geom_vline(xintercept = c(-0.5, 0.5), linetype = 2) +
  geom_hline(yintercept = -log10(0.05), linetype = 2) +
  coord_cartesian(xlim = x_lim, ylim = y_lim) +
  facet_grid(contrast ~ lfc_type, switch = "y") +
  scale_colour_manual(values = c("FALSE" = "grey60", "TRUE" = "steelblue4"),
                      labels = c("FALSE" = "FDR>0.05 or LFC<0.5", "TRUE" = "FDR<=0.05 &  LFC>=0.5"),name = NULL) +
  labs(
    title = "Volcano comparison: Raw vs Shrunken (ashr) with LFC filtering",
    x = "log2 fold change",
    y = expression(-log[10]("FDR from unshrunken"))) +
  theme_minimal(base_size = 11) +
  theme(
    legend.position = "bottom",
    strip.placement = "outside",
    strip.background = element_blank(),
    strip.text.x = element_text(face = "bold"),
    strip.text.y.left = element_text(angle = 0, face = "bold"),
    panel.grid.minor = element_blank())

## Without LFC filtering

# Raw vs Shrunken
res_long_nolfc<- res_all %>%
  mutate(
    sig_raw = !is.na(FDR_raw) & FDR_raw <= 0.05,
    sig_shr = !is.na(FDR_raw) & FDR_raw <= 0.05
  ) %>%
  dplyr::select(gene, contrast, FDR_raw, lfc_raw, lfc_shr, sig_raw, sig_shr) %>%
  tidyr::pivot_longer(
    cols= c(lfc_raw, lfc_shr, sig_raw, sig_shr),
    names_to= c(".value", "lfc_type"),
    names_sep= "_") %>%
  mutate(
    lfc_type = recode(lfc_type,
      raw = "Raw LFC (before)",
      shr = "Shrunken LFC (after)"),
    lfc_type = factor(lfc_type, levels = c("Raw LFC (before)", "Shrunken LFC (after)")),
    contrast= factor(contrast, levels = c(
      "Female: Ang-II vs Saline",
      "Male: Ang-II vs Saline",
      "Interaction: (M & F) treatment effect",
      "Baseline (Saline): M vs F"
    )))

x_max<- quantile(abs(res_long_nolfc$lfc), 0.995, na.rm = TRUE)
x_lim<- c(-max(1.5, x_max), max(1.5, x_max))
y_max<- quantile(-log10(res_long_nolfc$FDR_raw), 0.995, na.rm = TRUE)
y_lim<- c(0, max(5, y_max))

ggplot(res_long_nolfc, aes(lfc, -log10(FDR_raw), colour = sig)) +
  geom_point(alpha = 0.6, size = 1.2) +
  geom_hline(yintercept = -log10(0.05), linetype = 2) +
  coord_cartesian(xlim = x_lim, ylim = y_lim) +
  facet_grid(contrast ~ lfc_type, switch = "y") +
  scale_colour_manual(
    values = c("FALSE" = "grey60", "TRUE" = "steelblue4"),
    labels = c("FALSE" = "FDR > 0.05", "TRUE" = "FDR ≤ 0.05"),
    name = NULL) +
  labs(
    title = "Volcano comparison: Raw vs Shrunken (ashr) without LFC filtering",
    x = "log2 fold change",
    y = expression(-log[10]("FDR (from unshrunken fit)"))) +
  theme_minimal(base_size = 11) +
  theme(
    legend.position = "bottom",
    strip.placement = "outside",
    strip.background = element_blank(),
    strip.text.x = element_text(face = "bold"),
    strip.text.y.left = element_text(angle = 0, face = "bold"),
    panel.grid.minor = element_blank())

Diagnostics

The interaction coefficient (SexM.TreatmentAng.II) is a difference-in-differences. If β_int ~ 0 for most genes, it means the Ang-II effect is about the same in males and females. Therefore, it could be assume that there can be still many DE genes within each sex but with small effects, however, the extra difference between males and females could almost be negligible for Ang-II.

Reduced model The LRT is comparing the full model to the reduced model to identify significant genes. The p-values are determined solely by the difference in deviance between the ‘full’ and ‘reduced’ model formula (not log2 fold changes). Essentially the LRT test is testing whether the term(s) removed in the ‘reduced’ model explains a significant amount of variation in the data. Quoted https://hbctraining.github.io/DGE_workshop/lessons/08_DGE_LRT.html

Context:

So in the case, where females are baseline/ reference; * β_int > 0: Ang-II effect is larger in males than females.

β_int < 0: Ang-II effect is smaller in males (or larger in females).

Results of the preliminary states that;

Group sizes: F (4 Ang-II, 3 Saline), M (4 Ang-II, 6 Saline) = unbalanced and small
F vs M LFC scatter: global Pearson r ≈ 0.056 is small but could be because ~95% of genes are near 0, hence noise dominates the correlation.
Summary of the full model (with interaction) to the reduced model (without interaction) gene-by-gene, shows ~0–4 genes at FDR <0.1 (and 0–1 at FDR < 0.05) out of 41594 genes. Hence it could be deduced that the interaction term adds almost no explanatory power overall.
Effect-size difference (M & F) distribution of median −0.023, centred near 0, further proving that Ang-II effect size is similar in males and females overall.
The histogram of LFC shows that the values of the res_F_raw and res_M_raw is centred around zero(median −0.023), indicating that Ang-II effects are broadly similar between sexes, consistent with the likelihood-ratio test showing negligible evidence for an interaction

# No.of saline and treatment samples of each sex
table(pheno$Sex, pheno$Intervention)

##    
##     Ang-II Saline
##   F      4      3
##   M      4      6

# correlation of the interaction 
summary(res_F_raw$lfcSE)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.05636 0.17710 0.40386 0.67548 0.90210 4.80267

summary(res_M_raw$lfcSE)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.04774 0.15084 0.35101 0.58459 0.78546 4.05917

cor.test(res_F_raw$log2FoldChange, res_M_raw$log2FoldChange, use="complete.obs")

## 
##  Pearson's product-moment correlation
## 
## data:  res_F_raw$log2FoldChange and res_M_raw$log2FoldChange
## t = 11.464, df = 41592, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.04653857 0.06569855
## sample estimates:
##        cor 
## 0.05612373

lfc_F<- res_F_raw$log2FoldChange         
lfc_M<- res_M_raw$log2FoldChange       
delta<- lfc_M - lfc_F   
summary(delta)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## -30.00000  -0.32066  -0.02349  -0.03099   0.33085  21.20516

# Model without interaction term (Hypothesis testing: Likelihood ratio test (LRT))
dds_lrt<- DESeq(dds, test = "LRT", reduced = ~ Sex + Treatment)
lrt<- results(dds_lrt)   
summary(lrt)

## 
## out of 41594 with nonzero total read count
## adjusted p-value < 0.1
## LFC > 0 (up)       : 3, 0.0072%
## LFC < 0 (down)     : 1, 0.0024%
## outliers [1]       : 21, 0.05%
## low counts [2]     : 0, 0%
## (mean count < 0)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results

lrt05<- results(dds_lrt, alpha = 0.05)
summary(lrt05)

## 
## out of 41594 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up)       : 0, 0%
## LFC < 0 (down)     : 1, 0.0024%
## outliers [1]       : 21, 0.05%
## low counts [2]     : 0, 0%
## (mean count < 0)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results

# Distribution of M & F interaction LFC
hist(res_int_raw$log2FoldChange, breaks=80)

Heatmaps of top 40 genes per contrast

Interaction heatmap shows the “top 40 by smallest FDR”

top_heatmap<- function(res_shr, vsd, label, n=40) {
  keep<- order(res_shr$padj, na.last = NA)[1:min(n, sum(!is.na(res_shr$padj)))]
  genes<- rownames(res_shr)[keep]
  mat<- assay(vsd)[genes, , drop=FALSE]
  mat<- scale(t(scale(t(mat))))  
  ann<- as.data.frame(colData(vsd))[, c("Sex","Treatment"), drop=FALSE]
  pheatmap(mat, annotation_col = ann, show_rownames = TRUE,
           main = paste0("Top ", n, " DE genes: ", label))
}

top_heatmap(res_F_shr, vsd, "F Ang-II vs Saline")

top_heatmap(res_M_shr, vsd, "M Ang-II vs Saline")

top_heatmap(res_int_shr, vsd, "Interaction (M & F)")

top_heatmap(res_bsl_shr, vsd, "Baseline (M vs F, Saline)")

# Female-only view of female contrast genes
female_cols<- colData(vsd)$Sex == "F"
top_heatmap(res_F_shr, vsd[, female_cols], "F Ang-II vs Saline (female samples)")

# Male-only view of male contrast genes
male_cols<- colData(vsd)$Sex == "M"
top_heatmap(res_M_shr, vsd[, male_cols], "M Ang-II vs Saline (male samples)")

top_heatmap_sig<- function(res_shr, vsd, label, n=40, alpha=0.05){
  sig<- which(!is.na(res_shr$padj) & res_shr$padj <= alpha)
  if (length(sig) == 0) {
    message("No genes pass FDR <= ", alpha, " for: ", label)
    return(invisible(NULL))
  }
  keep<- head(sig[order(res_shr$padj[sig])], n)
  genes<- rownames(res_shr)[keep]
  mat<- assay(vsd)[genes, , drop=FALSE]
  mat<- scale(t(scale(t(mat))))
  ann<- as.data.frame(colData(vsd))[, c("Sex","Treatment"), drop=FALSE]
  pheatmap(mat, annotation_col=ann, show_rownames=TRUE,
           main=sprintf("Top %d genes (FDR <= %.2g): %s", length(genes), alpha, label))
}

top_heatmap_sig(res_F_shr, vsd, "F Ang-II vs Saline")

top_heatmap_sig(res_M_shr, vsd, "M Ang-II vs Saline")

top_heatmap_sig(res_bsl_shr, vsd, "Baseline (M vs F, Saline)")

sum(res_int_shr$padj <= 0.05, na.rm=TRUE)

## [1] 1

# Female-only view of female contrast genes
female_cols<- colData(vsd)$Sex == "F"
top_heatmap_sig(res_F_shr, vsd[, female_cols], "F Ang-II vs Saline (female samples)")

# Male-only view of male contrast genes
male_cols<- colData(vsd)$Sex == "M"
top_heatmap_sig(res_M_shr, vsd[, male_cols], "M Ang-II vs Saline (male samples)")

Pathway enrichment Analysis

GSEA (Hallmarks)

MSigDB “Hallmarks” (collection H) = 50 core biological with minimal overlap gene sets that capture broad, well-defined biological processes.

# Map Ensembl to gene symbols 
gtf<- file.path(idx, "Mus_musculus.GRCm39.115.gtf")
gtf_df<- as.data.frame(rtracklayer::import(gtf))

# For symbols
ann<- gtf_df %>%
  filter(type == "gene") %>%
  transmute(
    ensgene = sub("\\.\\d+$", "", gene_id),
    symbol  = gene_name) %>%
  distinct()

head(ann)

##              ensgene  symbol
## 1 ENSMUSG00000104478 Gm38212
## 2 ENSMUSG00000104385  Gm7449
## 3 ENSMUSG00000101231 Gm28283
## 4 ENSMUSG00000101097  Gm6679
## 5 ENSMUSG00000100764 Gm29155
## 6 ENSMUSG00000100831 Gm17847

# map ensembl to entrez
map_ensg_to_entrez<- function(ens_ids){
  ens_ids <- sub("\\.\\d+$","", ens_ids)
  tibble(ENSEMBL = ens_ids) %>%
    mutate(ENTREZID = AnnotationDbi::mapIds(
      org.Mm.eg.db, keys = ENSEMBL, keytype = "ENSEMBL",
      column = "ENTREZID", multiVals = "first")) %>%
    filter(!is.na(ENTREZID)) %>% distinct()
}

# rank vectors on entrez 
rankify_entrez<- function(res_raw, ann_map){
  df<- as.data.frame(res_raw) %>%
    rownames_to_column("ensgene") %>%
    mutate(ensgene = sub("\\.\\d+$","", ensgene)) %>%
    left_join(ann_map, by = c("ensgene"="ENSEMBL")) %>%
    filter(!is.na(stat), !is.na(ENTREZID))
  split(df$stat, df$ENTREZID) %>% vapply(\(v) v[which.max(abs(v))], numeric(1))
}

ens_all<- unique(c(rownames(res_F_raw), rownames(res_M_raw),
                       rownames(res_int_raw), rownames(res_bsl_raw)))
ann_entrez<- map_ensg_to_entrez(ens_all)

ranks_F_e<- rankify_entrez(res_F_raw, ann_entrez)
ranks_M_e<- rankify_entrez(res_M_raw, ann_entrez)
ranks_Int_e<- rankify_entrez(res_int_raw, ann_entrez)
ranks_Bas_e<- rankify_entrez(res_bsl_raw, ann_entrez)

# Hallmark pathways grouped by gs_name
msig_entrez<- msigdbr(species="Mus musculus", category="H") %>%
  distinct(gs_name, entrez_gene) %>%
  group_by(gs_name) %>%
  summarise(genes = list(unique(entrez_gene)), .groups = "drop") %>%
  tibble::deframe()

# overlap 
ov_Fe<- sapply(msig_entrez, function(g) sum(names(ranks_F_e) %in% g)); summary(ov_Fe)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   32.00   95.75  164.50  142.28  195.75  202.00

# GSEA
run_gsea<- function(ranks, label, pathways) {
  fgsea::fgseaMultilevel(pathways = pathways, stats = ranks, minSize = 10, maxSize = 2000) %>% arrange(padj) %>% mutate(contrast = label)
}

gsea_F<- run_gsea(ranks_F_e, "F Ang-II vs Saline", msig_entrez)
gsea_M<- run_gsea(ranks_M_e, "M Ang-II vs Saline", msig_entrez)
gsea_Int<- run_gsea(ranks_Int_e, "Interaction (M & F)", msig_entrez)
gsea_Bas<- run_gsea(ranks_Bas_e, "Baseline (M vs F)", msig_entrez)

gsea_all<- bind_rows(gsea_F, gsea_M, gsea_Int, gsea_Bas)

length(ranks_F_e)

## [1] 22219

sum(duplicated(names(ranks_F_e)))

## [1] 0

summary(ranks_F_e)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -5.4778 -0.5577  0.2257  0.2475  1.0306  6.2055